home *** CD-ROM | disk | FTP | other *** search
- *-------------------------------------------------------------------------------
- *-- Program...: STATS.PRG
- *-- Programmer: Ken Mayer (CIS: 71333,1030) and Jay Parsons (CIS: 70160,340)
- *-- Date......: 06/25/1992
- *-- Notes.....: Statistical Functions -- see README.TXT to include this
- *-- library file in your system.
- *-------------------------------------------------------------------------------
-
- FUNCTION Samplevar
- *-------------------------------------------------------------------------------
- *-- Programmer..: Jay Parsons (CIS: 70160,340)
- *-- Date........: 4/13/1992
- *-- Notes.......: Finds sample variance of specified field of the current
- *-- : database, using CALCULATE command.
- *-- : The CALCULATE command calculates the population variance,
- *-- : which is smaller by a factor of (n-1)/n.
- *-- :
- *-- Written for.: dBASE IV Version 1.5
- *-- Rev. History: Original function 1990.
- *-- : Modified to take optional parameter, 4/13/1992
- *-- Calls : None
- *-- Called by...: Any
- *-- Usage.......: Samplevar( <cField> [, <cClause> ] )
- *-- Example.....: ? Samplevar( "Balance", ".FOR..NOT. isblank( Balance )" )
- *-- Returns : a numeric or float value, the sample variance, or .F. if
- *-- : it cannot be calculated.
- *-- : If any of the numeric items are floats, the result will be.
- *-- Parameters..: cField, name of a numeric field of the current database
- *-- : for which to calculate the sample variance
- *-- : cClause, optional, a FOR, WHILE, TO, etc. clause
- *-------------------------------------------------------------------------------
- PARAMETERS cField, cCondition
- PRIVATE fVar, nCount, cCond
- IF pcount() = 2
- cCond = " "+ cCondition
- ELSE
- cCond = ""
- ENDIF
- CALCULATE VAR( &cField ), CNT() TO fVar, nCount &cCond
-
- RETURN iif( nCount > 1, fVar * nCount / ( nCount - 1 ), .F. )
- *-- Eof: Samplevar()
-
- FUNCTION Stny
- *-------------------------------------------------------------------------------
- *-- Programmer..: Jay Parsons (CIS: 70160,340)
- *-- Date........: 11/13/1990
- *-- Notes.......: Returns value of the standard normal distribution function
- *-- : given a number of standard deviations from the mean.
- *-- : This function is not useful alone. The standard normal
- *-- : distribution function is the familiar bell-shaped curve
- *-- : scaled so its mean is at 0, each standard deviation is 1
- *-- : and the total area under the curve is 1. The function
- *-- : Stnarea calls on this function to calculate the approximate
- *-- : area (a fraction equal to percent of the total) under the
- *-- : part of the curve lying betwen the mean and the given
- *-- : number of standard deviations.
- *-- :
- *-- Written for.: dBASE IV
- *-- Rev. History: None
- *-- Calls : None
- *-- Called by...: Any
- *-- Usage.......: Stny( <nDevs> )
- *-- Example.....: ? Stny( 1 )
- *-- Returns : numeric value of the function.
- *-- Parameters..: nDevs, standard deviations from the mean
- *-------------------------------------------------------------------------------
- PARAMETERS nDevs
-
- RETURN exp( -nDevs * nDevs / 2 ) / sqrt( 2 * pi() )
- *-- EoF: Stny()
-
- FUNCTION Stnarea
- *-------------------------------------------------------------------------------
- *-- Programmer..: Jay Parsons (CIS: 70160,340)
- *-- Date........: 11/13/1990
- *-- Notes.......: Area of the standard normal distribution function between
- *-- : mean and given number of standard deviations from the mean.
- *-- :
- *-- : What's it about? Well, College Board scores (originally)
- *-- : were based on a normal distribution with a mean of 500 and
- *-- : 100 points per standard deviation. Knowing that a 650
- *-- : score is 1.5 standard deviations from the 500 mean, we
- *-- : can calculate Stnarea( 1.5 ) as .4332. This tells us that
- *-- : 43.32% of the scores lie between 650 and 500. Since 50%
- *-- : lie below 500, a score of 650 beats 93.32% of the scores.
- *-- :
- *-- : The polynomial approximation used by this function is said
- *-- : to be accurate to .00001, 1/1000 of one percent. Remember
- *-- : to SET DECIMALS appropriately to view results.
- *-- :
- *-- Written for.: dBASE IV
- *-- Rev. History: None
- *-- Calls : Stny() Function in STATS.PRG
- *-- Called by...: Any
- *-- Usage.......: Stnarea( <nDevs> )
- *-- Example.....: ? Stnarea( 1.5 )
- *-- Returns : % of area between deviations given and the mean, 0<=a<.5.
- *-- Parameters..: nDevs, standard deviations from the mean
- *-------------------------------------------------------------------------------
- PARAMETERS nDevs
- PRIVATE nX, nV
- nX = abs( nDevs )
- nV = 1 / ( 1 + .33267 * nX )
-
- RETURN .5 - Stny( nX ) * ( .4361836 * nV - .1201676 * nV * nV ;
- + .937298 * nV * nV * nV )
- *-- EoF: Stnarea()
-
- FUNCTION Stnz
- *-------------------------------------------------------------------------------
- *-- Programmer..: Jay Parsons (CIS: 70160,340)
- *-- Date........: 11/13/1990
- *-- Notes.......: A lookup table to find the values of "z", standard
- *-- : deviations, corresponding to the most common areas inside a
- *-- : given number of tails of the normal distribution function.
- *-- :
- *-- : Used in testing confidence intervals. If a sample of
- *-- : light bulbs from a shipment shows an average life of 1150
- *-- : hours, and the criterion for rejection of the shipment is
- *-- : 95% confidence that the average life of all bulbs is less
- *-- : than (a single tail) 1200 hours, the value 1.64485 returned
- *-- : by this function is necessary to determine whether to
- *-- : reject the shipment or not.
- *-- :
- *-- : Values of "z" that are not found in the table can be found
- *-- : using Stndevs, below, but it is slow.
- *-- :
- *-- Written for.: dBASE IV
- *-- Rev. History: None
- *-- Calls : None
- *-- Called by...: Any
- *-- Usage.......: Stnz( <nProb>, <nTails> )
- *-- Example.....: ? Stnz( .95, 1 )
- *-- Returns : z, number of standard deviations from mean inside which
- *-- : ( or to the side of which includes the mean, if one tail)
- *-- : the given percentage of area will fall.
- *-- : Returns -1 if no entry in table.
- *-- Parameters..: nConf, confidence desired, 0 < nConf < 1
- *-- : nTails, 1 or 2 = number of tails of curve of interest
- *-------------------------------------------------------------------------------
- PARAMETERS nConf, nTails
- IF nTails # 1 .AND. nTails # 2
- RETURN -1
- ENDIF
- DO CASE
- CASE nConf = .95
- RETURN iif( nTails = 1, 1.64485, 1.96010 )
- CASE nConf = .99
- RETURN iif( nTails = 1, 2.32676, 2.57648 )
- CASE nConf = .995
- RETURN iif( nTails = 1, 2.57648, 2.80794 )
- CASE nConf = .999
- RETURN iif( nTails = 1, 3.09147, 3.29202 )
- OTHERWISE
- RETURN -1
- ENDCASE
-
- *-- EoF: Stnz()
-
- FUNCTION Stndiff
- *-------------------------------------------------------------------------------
- *-- Programmer..: Jay Parsons (CIS: 70160,340)
- *-- Date........: 4/13/1992
- *-- Notes.......: Determines whether hypothesis that sample of a given mean
- *-- is different from expected mean is justified.
- *--
- *-- If nPopstd, the standard deviation of the population, is
- *-- not known and nSample, the sample size, is greater than
- *-- 30, the sample standard deviation may be used for nPopstd.
- *--
- *-- This function assumes the population is large relative to
- *-- the sample or that the sampling is with replacement. If
- *-- neither is true, the right side of the expression in the
- *-- later return line should be multiplied by:
- *-- sqrt( ( nPop - nSample ) / ( nPop - 1 ) )
- *-- where nPop is the size of the population.
- *--
- *-- Do not use this with small samples, less than 20, because
- *-- the standard normal distribution is not sufficiently
- *-- accurate as an approximation of the distribution of sample
- *-- means in such a case. See "Student's T-distribution" in a
- *-- statistics text.
- *--
- *-- Written for.: dBASE IV Version 1.5
- *-- Rev. History: None.
- *-- Calls : Stnz() Function in STATS.PRG
- *-- Called by...: Any
- *-- Usage.......: Stndiff( <nConf>, <nTails>, <nSample>, <nSampmean>, ;
- *-- : <nPopmean>, <nPopstd> )
- *-- Example.....: ? Stndiff( .95, 1, 30, 1150, 1200, 20 )
- *-- Returns : .T. if hypothesis of difference is justified to degree of
- *-- : confidence specified, or .F. Returns -1 if confidence is
- *-- : not one for which z can be looked up in Stnz(). If you
- *-- : need other confidence levels, run Stndevs() to find the
- *-- : z values for them and add them to the Stnz() table.
- *-- Parameters..: nConf, confidence desired, 0 < nConf < 1
- *-- : nTails, 1 or 2 = number of tails of curve of interest
- *-- : nSample, number of items in the sample
- *-- : nSampmean, mean of the sample
- *-- : nPopmean, mean of the population ( test standard mean )
- *-- : nPopstd, standard deviation of population
- *-------------------------------------------------------------------------------
- PARAMETERS nConf, nTails, nSample, nSampmean, ;
- nPopmean, nPopstd
- PRIVATE nStd
- nStd = Stnz( nConf, nTails )
- IF nStd = -1
- RETURN nStd
- ELSE
- RETURN abs( nSampmean - nPopmean ) ;
- > nStd * nPopstd / sqrt( nSample )
- ENDIF
- *-- EoF: Stndiff()
-
- FUNCTION Stndevs
- *-------------------------------------------------------------------------------
- *-- Programmer..: Jay Parsons (CIS: 70160,340)
- *-- Date........: 4/13/1992
- *-- Notes.......: Calculates "z", standard deviations, corresponding to any
- *-- : area of standard normal curve between mean and the desired
- *-- : z. Much slower than Stnz().
- *-- Written for.: dBASE IV Version 1.5
- *-- Rev. History: Original function 1990.
- *-- : Conformed to Zeroin() 4/13/1992.
- *-- Calls : Zeroin() Function in STATS.PRG
- *-- Called by...: Any
- *-- Usage.......: Stndevs( <nArea> )
- *-- Example.....: ? Stndevs( .96 )
- *-- Returns : z, number of standard deviations from mean, or a negative
- *-- : number indicating failure to find a root..
- *-- Parameters..: nArea, area "left" of point of interest, .5 < nArea < 1
- *-------------------------------------------------------------------------------
- PARAMETERS nArea
- PRIVATE nTest, nFlag
- IF nArea > .99999 .OR. nArea < .5
- RETURN -1
- ENDIF
- nFlag = 0
- nTest = Zeroin( "Tstnarea", 0, 4.2, float(1/100000), 100, nFlag, nArea )
-
- RETURN iif( nFlag < 3, nTest, -nFlag )
- *-- EoF: Stndevs()
-
- FUNCTION Tstnarea
- *-------------------------------------------------------------------------------
- *-- Programmer..: Jay Parsons (CIS: 70160,340)
- *-- Date........: 11/13/1990
- *-- Notes.......: Translation function to convert area to left of point
- *-- : under standard normal curve to 0 for Zeroin().
- *-- Written for.: dBASE IV
- *-- Rev. History: None
- *-- Calls : Stnarea() Function in STATS.PRG
- *-- Called by...: Any
- *-- Usage.......: Tstnarea( <nDevs>, <nArea> )
- *-- Example.....: ? Tstnarea( 1.6,.96 )
- *-- Returns : positive or negative number corresponding to direction to
- *-- : root where nArea = Stnarea
- *-- Parameters..: nDevs, trial number of standard deviations
- *-- : nArea, area for which deviations are to be found
- *-------------------------------------------------------------------------------
- PARAMETERS nDevs, nArea
-
- RETURN Stnarea( nDevs ) +.5 - nArea
- *-- EoF: Tstnarea()
-
- FUNCTION Zeroin
- *-------------------------------------------------------------------------------
- *-- Programmer..: Tony Lima (CIS: 72331,3724) and Jay Parsons (CIS: 70160,340)
- *-- Date........: 04/13/1992
- *-- Notes.......: Finds a zero of a continuous function.
- *-- : In substance, what this function does is close in on a
- *-- : solution to a function that cannot otherwise be solved.
- *-- : Assuming Y = f(X), if Y1 and Y2, the values of the function
- *-- : for X1 and X2, have different signs, there must be at least
- *-- : one value of X between X1 and X2 for which Y = 0, if the
- *-- : function is continuous. This function closes in on such a
- *-- : value of X by a trial-and-error process.
- *-- :
- *-- : This function is very slow, so a maximum number of iterations
- *-- : is passed as a parameter. If the number of iterations is
- *-- : exceeded, the function will fail to find a root. If this
- *-- : occurs, pick different original "X" values, increase the
- *-- : number of iterations or increase the errors allowed. Once
- *-- : an approximate root is found, you can use values of X close
- *-- : on either side and reduce the error allowed to find an
- *-- : improved solution. Also, of course, the signs of Y must be
- *-- : different for the starting X values for the function to
- *-- : proceed at all.
- *-- :
- *-- : NOTE ESPECIALLY - There is NO guarantee that a root returned
- *-- : by this function is the only one, or the most meaningful.
- *-- : It depends on the function that this function calls, but if
- *-- : that function has several roots, any of them may be returned.
- *-- : This can easily happen with such called functions as net
- *-- : present value where the cash flows alternate from positive
- *-- : to negative and back, and in many other "real life" cases.
- *-- : See the discussion of @IRR in the documentation of a good
- *-- : spreadsheet program such as Quattro Pro for further
- *-- : information.
- *-- :
- *-- : The method used by this function is a "secant and bisect"
- *-- : search. The "secant" is the line connecting two X,Y
- *-- : points on a graph using standard Cartesian coordinates.
- *-- : Where the secant crosses the X axis is the best guess for
- *-- : the value of X that will have Y = 0, and will be correct
- *-- : if the function is linear between the two points. The
- *-- : basic strategy is to calculate Y at that value of X, then
- *-- : keep the new X and that one of the old X values that had
- *-- : a Y-value of opposite sign, and reiterate to close in.
- *-- :
- *-- : If the function is a simple curve with most of the change
- *-- : in Y close to one of the X-values, as often occurs if the
- *-- : initial values of X are poorly chosen, repeated secants
- *-- : will do little to find a Y-value close to zero and will
- *-- : reduce the difference in X-values only slightly. In this
- *-- : case the function shifts to choosing the new X halfway
- *-- : between the old ones, bisecting the difference and always
- *-- : reducing the bracket by half, for a while.
- *-- :
- *-- : While this function finds a "zero", it may be used to
- *-- : find an X corresponding to any other value of Y. Suppose
- *-- : the function of X is FUNCTION Blackbox( X ) and it is
- *-- : desired to find a value of X for which f(X) = 7. The trick
- *-- : is to interpose a function between Zeroin() and Blackbox()
- *-- : that will return a 0 to Zeroin() whenever Blackbox() returns
- *-- : 7. By calling that function, Zeroin() finds a value of
- *-- : X for which Blackbox( X ) = 7, as required:
- *-- : Result = Zeroin( "Temp", <other parameters omitted> )
- *-- :
- *-- : FUNCTION Temp
- *-- : parameters nQ
- *-- : RETURN Blackbox( nQ ) - 7
- *-- :
- *-- Written for.: dBASE IV Version 1.5
- *-- Rev. History: Original function 1990.
- *-- : Modified to take optional parameters, 4/13/1992
- *-- Calls : The function whose name is first parameter.
- *-- : NPV() Function in FINANCE.PRG
- *-- Called by...: Any
- *-- Usage.......: Zeroin( <cFunction>, <fX1>, <fX2>, <fAbserror>, ;
- *-- : <nMaxiter>, <n_Flag> ;
- *-- : [, xPass1 [, xPass2 [, xPass3 ] ] ] )
- *-- Example.....: ? Zeroin( "Npv", 0, 200, .000001, 200, n_Flag, 11 )
- *-- Returns : a float value representing a root, if n_Flag < 3.
- *-- Parameters..: cFunction, the name of the function to solve for a root.
- *-- fX1, one of the X-values between which the root is sought.
- *-- fX2, the second of these values.
- *-- Note: These MUST be chosen so the f( X ) values for the two
- *-- of them have opposite signs (they must bracket the result).
- *-- fAbserror, the absolute error allowed in the result.
- *-- nMaxiter, the maximum number of times to iterate.
- *-- n_Flag, an integer to signal success ( < 3 ) or failure.
- *-- xPass1 . . . 3, arguments to be passed through to cFunction.
- *-- The parameter "n_Flag" should be passed as a variable so it
- *-- may be accessed on return. The limit of 9 literal parameters
- *-- may require passing others as variables. The "xPass"
- *-- parameters are optional and the fact there are three of them
- *-- is arbitrary; they exist to hold whatever parameters may be
- *-- needed by the function cFunction being called aside from
- *-- the value of X for which it is being evaluated. Add more
- *-- and change the 3 "&cFunc." lines below if you need more.
- *-- Side effects: Uses and alters a global numeric variable, here called
- *-- "n_Flag", to report error conditions resulting in value
- *-- returned being meaningless. Possible n_Flag values are:
- *-- 1 success - root found within error allowed
- *-- 2 success - root was found exactly
- *-- 3 error - function value not converging
- *-- 4 error - original values do not bracket a root
- *-- 5 error - maximum iterations exceeded
- *-------------------------------------------------------------------------------
- parameters cFunc, fNearx, fFarx, fAbserr, nMaxiter, ;
- n_Flag, xPass1, xPass2, xPass3
- private nSplits, fBracket, fFary, fNeary, nIters
- private fMaxabs, fOldx, fOldy, fDiffx, fAbsdiff, fSecant
-
- store 0 to nSplits, nIters
- fBracket = abs ( fNearx - fFarx )
- fFary = &cFunc.( fFarx, xPass1, xPass2, xPass3 )
- fNeary = &cFunc.( fNearx, xPass1, xPass2, xPass3 )
-
- if sign( fNeary ) = sign( fFary )
- n_Flag = 4
- return float(0)
- endif
-
- fMaxabs = max( abs( fNeary ), abs( fFary ) )
- n_Flag = 0
-
- * Main iteration loop
-
- do while .t.
-
- if abs( fFary ) < abs( fNeary )
-
- * Interchange fNearx and fFarx so that
- * fNearx is closer to a solution--
- * abs( fNeary ) <= abs( fFary )
-
- fOldx = fNearx
- fOldy = fNeary
- fNearx = fFarx
- fNeary = fFary
- fFarx = fOldx
- fFary = fOldy
- endif
-
- fDiffx = fFarx - fNearx
- fAbsdiff = abs( fDiffx )
-
- * Test whether interval is too small to continue
-
- if fAbsdiff <= 2 * fAbserr
- if abs( fNeary ) > fMaxabs
-
- * Yes, but we are out of bounds
-
- n_Flag = 3
- fNearx = float(0)
- else
-
- * Yes, and we have a solution!
-
- n_Flag = 1
- endif
- exit
- endif
-
- * Save the last approximation to x and y
-
- fOldx = fNearx
- fOldy = fNeary
-
- * Check if reduction in the size of
- * bracketing interval is satisfactory.
- * If not, bisect until it is.
-
- nSplits = nSplits + 1
- if nSplits >= 4
- if 4 * fAbsdiff >= fBracket
- fNearx = fNearx + fDiffx / 2
- else
- nSplits = 0
- fBracket = fAbsdiff / 2
-
- * Calculate secant
-
- fSecant = ( fNearx - fFarx ) * fNeary ;
- / ( fFary - fNeary )
-
- * But not less than error allowed
-
- if abs( fSecant ) < fAbserr
- fNearx = fnearx + fAbserr * sign( fDiffx )
- else
- fNearx = fNearx + fSecant
- endif
- endif
- endif
-
- * Evaluate the function at the new approximation
-
- fNeary = &cFunc.( fNearx, xPass1, xPass2, xPass3 )
-
- * If it's exactly zero, we win! Run with it
-
- if fNeary = 0.00
- n_Flag = 2
- exit
- endif
-
- * Else adjust iteration count and quit if too
- * many iterations with no solution
-
- nIters = nIters + 1
- if nIters > nMaxiter
- n_Flag = 5
- fNearx = float( 0 )
- exit
- endif
-
- * And finally keep as the new fFarx that one
- * of the previous approximations, fFarx and
- * fOldx, at which the function has a sign opposite
- * to that at the new approximation, fNearx.
-
- if sign( fNeary ) = sign( fFary )
- fFarx = fOldx
- fFary = fOldy
- endif
- enddo
-
- RETURN fNearx
- *-- EoF: Zeroin()
-
- *-------------------------------------------------------------------------------
- *-- The functions below are here by courtesy ... (to make life easier on the
- *-- poor programmer ...)
- *-------------------------------------------------------------------------------
-
- FUNCTION Npv
- *-------------------------------------------------------------------------------
- *-- Programmer..: Tony Lima (CIS: 72331,3724) and Jay Parsons (CIS: 70160,340)
- *-- Date........: 03/01/1992
- *-- Notes.......: Net present value of array aCashflow[ nPeriods ]
- *-- Calculates npv given assumed rate and # periods.
- *-- Written for.: dBASE IV, 1.1
- *-- Rev. History: None
- *-- Calls.......: None
- *-- Called by...: Any
- *-- Usage.......: NPV(<nRate>,<nPeriods>)
- *-- Example.....: ? NPV( .06, 6 )
- *-- Returns.....: Float = value of the project at given rate
- *-- Parameters..: nRate = Interest Rate
- *-- : nPeriods = Number of Periods to calculate for
- *-- Other inputs: Requires the array aCashflow[ ] set up before calling.
- *-- : Each of its elements [n] holds the cash flow at the
- *-- : beginning of period n, with a negative amount indicating
- *-- : a cash outflow. Elements of value 0 must be included for
- *-- : all periods with no cash flow, and all periods must be of
- *-- : equal length.
- *-- : If the project is expected to require an immediate outlay
- *-- : of $6,000 and to return $2,000 at the end of each of the
- *-- : first five years thereafter, the array will be:
- *-- : aCashflow[1] = -6000
- *-- : aCashflow[2] = 2000
- *-- : aCashflow[3] = 2000
- *-- : * * *
- *-- : aCashflow[6] = 2000
- *-- : Rewriting function to have array name passed as a parameter
- *-- : is possible, but will slow down execution to an extent that
- *-- : will be very noticeable if this function is being repeatedly
- *-- : executed, as by Zeroin() to find an Internal Rate of Return.
- *-------------------------------------------------------------------------------
-
- parameters nRate, nPeriods
- private nDiscount, nFactor, nPeriod, nNpv
- nPeriod = 1
- nNpv = aCashflow[ 1 ]
- nDiscount = float( 1 )
- nFactor = 1 / ( 1 + nRate )
- do while nPeriod < nPeriods
- nPeriod = nPeriod + 1
- nDiscount = nDiscount * nFactor
- nNpv = nNpv + aCashflow[ nPeriod ] * nDiscount
- enddo
-
- RETURN nNpv
- *-- EoF: Npv()
-
- FUNCTION ArrayRows
- *-------------------------------------------------------------------------------
- *-- Programmer..: Jay Parsons (CIS: 70160,340)
- *-- Date........: 03/01/1992
- *-- Notes.......: Number of Rows in an array
- *-- Written for.: dBASE IV, 1.1
- *-- Rev. History: None
- *-- Calls.......: None
- *-- Called by...: Any
- *-- Usage.......: ArrayRows("<aArray>")
- *-- Example.....: n = ArrayRows("aTest")
- *-- Returns.....: numeric
- *-- Parameters..: aArray = Name of array
- *-------------------------------------------------------------------------------
-
- parameters aArray
- private nHi, nLo, nTrial, nDims
- nLo = 1
- nHi = 1170
- if type( "&aArray[ 1, 1 ]" ) = "U"
- nDims = 1
- else
- nDims = 2
- endif
- do while .T.
- nTrial = int( ( nHi + nLo ) / 2 )
- if nHi < nLo
- exit
- endif
- if nDims = 1 .and. type( "&aArray[ nTrial ]" ) = "U" .or. ;
- nDims = 2 .and. type( "&aArray[ nTrial, 1 ]" ) = "U"
- nHi = nTrial - 1
- else
- nLo = nTrial + 1
- endif
- enddo
-
- RETURN nTrial
- *-- EoF: ArrayRows()
-
- *-------------------------------------------------------------------------------
- *-- End of Program: STATS.PRG
- *-------------------------------------------------------------------------------